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Abstract. - Many epidemic processes in networks spread by stochastic contacts among their 
connected vertices. There are two limiting cases widely analyzed in the physics literature, the 
so-called contact process (CP) where the contagion is expanded at a certain rate from an infected 
vertex to one neighbor at a time, and the reactive process (RP) in which an infected individual 
effectively contacts all its neighbors to expand the epidemics. However, a more realistic sce- 
nario is obtained from the interpolation between these two cases, considering a certain number of 
stochastic contacts per unit time. Here we propose a discrete-time formulation of the problem of 
contact-based epidemic spreading. We resolve a family of models, parameterized by the number 
of stochastic contact trials per unit time, that range from the CP to the RP. In contrast to the 
common heterogeneous mean-field approach, we focus on the probability of infection of individual 
nodes. Using this formulation, we can construct the whole phase diagram of the different infection 
models and determine their critical properties. 



The problem of modeling how diseases spread among 
individuals has been intensively studied for many years 
[HI]. The development of mathematical models to guide 
our understanding of the disease dynamics has allowed to 
address important issues such as immunization and vacci- 
nation policies [2J06]. Physicist's approaches to problems 
in epidemiology involve statistical physics, the theory of 
phase transitions and critical phenomena [7], which have 
been extremely helpful to grasp the macroscopic behavior 
of epidemic outbreaks [8TfT6] . The main artifice of this suc- 
cess has been the Mean-Field (MF) approximation, where 
local homogeneities of the ensemble are used to average 
the system, reducing degrees of freedom. 



The study of complex networks 1 7 1 !)] has provided 
new grounds to the understanding of contagion dynam- 
ics. Particularly important in nature are scale-free (SF) 
networks, whose degree distribution follows a power law 
P(k) ~ fc~ 7 for the number of connections, k, an individ- 



ual has. SF networks include patterns of sexual contacts 
[20] , the Internet [21] , as well as other social, technological 
and biological networks [35]. The critical properties of an 
epidemic outbreak in SF networks can be addressed using 
the heterogeneous MF (HMF) prescription [8-15]. It con- 
sists of coarse-grained vertices within degree classes and 
considers that all nodes in a degree class have the same 
dynamical properties; the approach also assumes that fluc- 
tuations can be neglected. Specifically, if j3 is the rate 
(probability per unit time) at which the disease spreads, 
it follows that the epidemic threshold in uncorrelated SF 
networks is given [8] by f3 c = (k)/(k 2 ), leading to (3 C — > 
as JV -> 00 when 2 < 7 < 3. 

MF approaches are extremely useful to assess the crit- 
ical properties of epidemic models, however they are not 
designed to give information about the probability of indi- 
vidual nodes but about classes of nodes. Then, questions 
concerning the probability that a given node be infected 
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are not well posed in this framework. To obtain more de- 
tails at the individual level of description, one has to rely 
on Monte Carlo (MC) simulations, which have also been 
used to validate the results obtained using MF methods. 
Restricting the scope of epidemiological models to those 
based in two states [313111] —susceptible (S) and infected 
(I) — , the current theory concentrates on two specific situ- 
ations, the contact process [23H28] (CP) and the reactive 
process [29 31J (RP). A CP stands for a dynamical pro- 
cess that involves an individual stochastic contagion per 
infected node per unit time, while in the RP there are as 
many stochastic contagions per unit time as neighbors a 
node has. This latter process underlies the abstraction 
of the susceptible-infected-susceptible (SIS) model [I][3j4]. 
However, in real situations, the number of stochastic con- 
tacts per unit time is surely a variable of the problem itself 

In this work, we introduce a theoretical framework for 
contact-based spreading of diseases in complex networks. 
Our formulation is based on probabilistic discrete-time 
Markov chains, generalizes existing HMF approaches and 
applies to weighted and unweighted complex networks. 
Within this context, in addition to capturing the global 
dynamics of the different contact models and its associated 
critical behavior, it is now possible to quantify the micro- 
scopic dynamics at the individual level by computing the 
probability that any node is infected in the asymptotic 
regime. MC simulations corroborate that the formalism 
introduced here reproduces correctly the whole phase di- 
agram for model and real-world networks. Moreover, we 
capitalize on the new approach to address how the spread- 
ing dynamics depends on the number of contacts actually 
used by a node to propagate the disease. 

Contact-based epidemic spreading models. 

Let us consider a network made up of N nodes, whose 
connections are represented by the entries {oij } of an N- 
by-N adjacency matrix A. Additionally, in the most gen- 
eral case in which the network is weighted, we denote by 
{coij} the weights of the connections between nodes, being 
Wi = J2j w ij the total strength [33) of node i. The above 
quantities completely define the structure of the underly- 
ing graph. As for the dynamics, we consider a discrete 
two-state (S and I) contact-based process. Each node of 
the network represents an individual (or a place, a city or 
airport for example) and each edge is a connection along 
which the infection spreads. At each time step, an infected 
node makes a number A of trials to transmit the disease 
to its neighbors with probability /3 per unit time. This 
forms a Markov chain where the probability of a node be- 
ing infected depends only on the last time step. After some 
transient time, the previous dynamics sets the system into 
a stationary state in which the average density of infected 
individuals, p, defines the prevalence of the disease. 

We next look at the probability that any given node i 
is infected at the stationary state. We denote by the 
probability that a node i is in contact with a node j, defin- 



ing a matrix R. These entries represent the probabilities 
that existing links in the network are used to transmit the 
infection. If i and j are not connected, then = 0. Be- 
sides, /i stands for the rate at which infected nodes are 
recovered and get back to the susceptible class; and fi- 
nally, Pi(t) is the probability that a node i is infected at 
time t. With these definitions, the discrete-time version of 
the evolution of the probability of infection of any node i 
reads 

Pi (t+1) = {1 - qi {t))(l - Pi (t)) + (1 - fi) Pi (t) 

+ H(l - qi(t))Pi(t) (1) 

where qi (t) is the probability of node i not being infected 
by any neighbor 

N 

%(t)=n( i -^ft(*)) ( 2 ) 

3=1 

The first term on the right hand side of eq. ([IJ is the 
probability that node i is susceptible (1 — Pi(t)) and is 
infected (1 — Qi{t)) by at least a neighbor. The second 
term stands for the probability that node i is infected at 
time t and does not recover, and finally the last term takes 
into account the probability that an infected node recovers 
(fipi(t)) but is re-infected by at least a neighbor (1 — (£)) . 
Within this formulation, we are assuming the most gen- 
eral situation in which recovery and infection occur on the 
same time scales, allowing then reinfection of individu- 
als during a discrete time window (for instance, one MC 
step). This formulation generalizes previous approxima- 
tions where reinfections can not occur. 

The formulation so far relies on the assumption that the 
probabilities of being infected pi are independent random 
variables. This hypothesis turns out to be valid in the vast 
majority of complex networks because the inherent topo- 
logical disorder makes dynamical correlations not persis- 
tent. The dynamical system ([l] [2]) corresponds to a family 
of possible models, parameterized by the explicit form of 
the contact probabilities r.y . Without loss of generality, it 
is instructive to think of these probabilities as the transi- 
tion probabilities of random walkers on the network. The 
general case is represented by Xi random walkers leaving 
node i at each time step: 



The CP corresponds to a model dynamics of one contact 
per unit time, Xi = 1, Vi in eq. (J3J thus ry = ufy/iujj. In 
the RP all neighbors are contacted, which corresponds, in 
this description, to set the limit A^ — » oo, Vz resulting on 
fij = aij regardless of whether the network is weighted or 

1 Strictly speaking, when A = 1, our model is not exactly the 
standard CP, since there reinfections are not considered. However, 
we will refer to it as a CP since only one neighbor is contacted at 
each time step and the critical points of both variants are the same. 
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not. Other prescriptions for Ai conform the spectrum of 
models that can be obtained using this unified framework. 
The phase diagram of every model is simply obtained solv- 
ing the system formed by eq. (HJ for i = 1, . . . , TV at the 
stationary state 

Pi = (1 - ®) + (1 - ( 4 ) 

This equation has always the trivial solution pi = 0, 
Vi = 1,...,N. Other non-trivial solutions are reflected 
as non zero fixed points of eq. ((4]) and can be easily com- 
puted numerically by iteration. The macroscopic order 
parameter is given by the expected infection density p, 
computed as 

1 N 

i=l 

Numerical results. — To show the validity of the ap- 
proach here discussed, we have performed MC simulations 
on different SF networks for RP. Figure [T] shows a com- 
parison of the phase diagram of the system obtained by 
MC simulations, with the numerical solution of eq. (j4]). To 
model the epidemic dynamics on the described topologies 
we incorporate a SIS model in which, at each time step, 
each node can be susceptible or infected. In our simula- 
tions time is discretized in time-steps and each simulation 
starts with a fraction po of randomly chosen infected indi- 
viduals (po = 0.05 in our simulations). At each time step 
an infected node i infects with the same probability /3 all 
its neighbors and recovers at a rate p. The simulation 
runs until a stationary state for the density of suscepti- 
ble individuals, p(t) is reached. The agreement between 
both curves is matchless. Moreover, the formalism also 
captures the microscopic dynamics as given by the p^s, 
see the inset of fig. [U 

In Figure 2 we analyze our formalism on top of the air- 
ports network data set, composed of passenger flights op- 
erating in the time period November 1, 2000, to October 
31, 2001 compiled by OAG Worldwide (Downers Grove, 
IL) and analyzed previously by Prof. Amaral's group [33] . 
It consists of 3618 nodes (airports) and 14142 links, we 
used the weighted network in our analysis. Airports cor- 
responding to a metropolitan area have been collapsed into 
one node in the original database. We show the density 
of infected individuals p as a function of j3 for different 
values of A. The critical points as well as the shape of 
the p — j3 phase diagrams greatly change at varying the 
number of stochastic contacts (A). For small values of A 
the disease prevalence is moderate, even for large values 
of the spreading rate /3. In contrast, when the number of 
trials is of order 10 3 the situation is akin to a RP. 

Finally, we compare the results of the formalism for dif- 
ferent random scale-free networks satisfying P(k) ~ /c~ 7 
generated using the configuration model. Figure |3] shows 
the phase diagram for p = 1 and several values of the 
exponent 7, both below and above 7 = 3. The system 
size has been fixed to N = 10 4 nodes. The dotted lines 




Fig. 1: Average fraction of infected individuals, p, as a function 
of the spreading rate /3 for N = 10 4 . The symbols correspond 
to MC simulations of the SIS model on top of random scale-free 
networks with 7 = 2.7 (error bars are smaller than the size of 
the symbol) and the lines stand for the analytical solutions of 
our formalism (with A = 00). We also represent in the inset a 
scatter plot for the probability that a node i (i — 1, . . . , N) is 
infected using results of MC simulations (the y-axis) and the 
solutions (x-axis) of eq. (|4]). Both results have been obtained 
for p = 1, the inset is for f3 — 0.1. 

represent the results obtained using the analytical approx- 
imation while symbols stand for MC simulations. As it can 
be seen, the agreement between both methods is remark- 
able, even for values of 7 < 2.5 where structural changes 
are extremely relevant 34 . The same agreement between 
MC results and the analytical solutions is obtained if one 
instead fixes the degree distribution exponent 7 and ex- 
plores the dependency with the system size. This is what 
is shown in fig. where we have depicted the phase di- 
agram for networks with 7 = 2.7 for several system sizes 
ranging from N = 500 to N = 10 5 . Except for N = 500, 
where MC results have a large standard deviation close to 
the critical point, the agreement is again excellent in the 
whole range of /3 values. 

Epidemic Threshold. — Let us now assume the ex- 
istence of a critical point /3 C for fixed values of p and Aj 
such that p = if (3 < (3 C and p > when (3 > (3 C . The 
calculation of this critical point is performed by consider- 
ing that when j3 — > /3 C , the probabilities pi rs ej, where 
< ti <C 1, and then after substitution in eq. ([2]) one gets 

N 

(6) 

i=i 

Inserting eq. © in eq. (U), and neglecting second order 
terms in e we get 




Vi = l,...,N (7) 
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Fig. 2: Density of infected individuals p as a function of /3 for 
different values of A in the air transportation network [33] . We 
have set p — 1 and p is calculated according to eq. ((5| once the 
Pi's are obtained. 

where stands for the Kronecker delta. The system ([7]) 
has non trivial solutions if and only if fj,//3 is an eigenvalue 
of the matrix R. Since we are looking for the onset of the 
epidemic, the lowest value of (3 satisfying ([7]) is 



(8) 



where A max is the largest eigenvalue of the matrix R. 
Equation ([5]) defines the epidemic threshold of the disease 
spreading process. 

It is worth analyzing the two limiting cases of CP and 
RP above. In the first case, one obtains the trivial result 
that the only non-zero solution corresponds to /3 C = /i, 
because the matrix R is a transition matrix whose max- 
imum eigenvalue is always A max = 1. For the RP corre- 
sponding to the SIS spreading process usually adopted [8] , 
the classical result for uncorrelated SF networks is recov- 
ered because, in this case, the largest eigenvalue [35036] is 
A max = (k 2 )/{k). 

Mesoscopic equations at the critical point. 

Once the general framework given by the dynamical sys- 
tem (fi] [2|) has been proposed, it is instructive to approxi- 
mate it using the hypotheses underlying HMF. These hy- 
potheses consist of: i) coarse-graining the system in classes 
of node by degree, assuming that the dynamical properties 
within each class are the same, and ii) neglecting fluctu- 
ations. To obtain the mesoscopic description we consider 
the second order approximation of eq. (j4|) and proceed as 
in the previous section. Therefore, 

qt « 1 - rjiCj + /3 2 rjirucjet (9) 

3 3<l 

After substitution in (QJ and reordering terms one gets 



= 



-fiei + /3(1 - £i) y rjiej + /Lt/3e 



i 
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Fig. 3: Phase diagram for the SIS model (A = oo) in a random 
scale free network for different 7's. The networks are made up 
oi N = 10 4 nodes and fx = 1. MC results are averages over 
10 2 realizations. Dashed lines corresponds to the theoretical 
prediction and symbols to MC results. 
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which are the equations governing the dynamics of the 
contact-based epidemic spreading process at the micro- 
scopic level. It is possible to write eq. (| 10[) at the com- 
monly used mesoscopic (degree class) level for unweighted, 
undirected heterogeneous networks. The interactions then 
takes place between classes of nodes. Denning the av- 
erage density of infected nodes with degree k as pk = 
j^J2k =kPi> wn ere Nh is the number of nodes with de- 
gree k and the sum runs over the set of nodes of degree k, 
we obtain the generalized HMF equation near criticality. 

Homogeneous networks. For homogeneous un- 
weighted undirected networks, ej = e and fc; « (k) for 
all nodes. Thus, p = Ylj e j — 6 an( l 



(3 2 p 2 ^2r Jiru 

3<l 



Defining 



(11) 
(12) 



R x (x) = 1 - (1 - x) x 
the terms involving values of rji are 

ra « ajiRxdk)- 1 ) (13) 



(k)Rx(ik)- 1 ) 



(14) 



J2 r ^ » ^(km-VRxttk)- 1 ) 2 (15) 

3<l 

Now, eq. ([TT]) becomes 

= -pp + Pp(l-p)(k)R x ((k)- v ) 
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Fig. 4: Phase diagram for the SIS model (A = oo) in a random 
scale free network for different system sizes as indicated. The 
networks have a power-law degree distribution with an expo- 
nent 7 = 2.7 and p, = 1. MC results are averages over 10 2 
realizations. 
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which may be considered as the MF approximation of our 
model for homogeneous networks. 

If A = 1 then Ri((k)~ l ) = -tK and eq. (fTrJf becomes 



(k) - 1 

Q=- f ,p + f3p(l-p) + f ,f3p 2 -^j^!3 2 p 2 (17) 

If A — > oo then i?oo((^} _1 ) = 1 and eq. (fl"6|) reads 

= -pp+Pp{l-p)(k)+p{3p 2 (k)- l -p 2 p 2 (k)({k)-l) (18) 

In both cases, the first two terms correspond to the 
standard CP and RP models (previously reported in the 
literature) respectively, and the additional terms are sec- 
ond order contributions corresponding to reinfections and 
multiple infections. 

Heterogeneous networks. Now we will concentrate on 
the class of heterogeneous unweighted undirected net- 
works completely specified by their degree distribution 
P(k) and by the conditional probability P(k'\k) that a 
node of degree k is connected to a node of degree k' . Of 
course, the normalization conditions P(k) = 1 and 
J2k> P{k'\k) — 1 must be fulfilled. In this case, the aver- 
age number of links that goes from a node of degree k to 
nodes of degree k' is kP(k'\k). 

In these heterogeneous networks it is supposed that all 
nodes of the same degree behave equally, thus £j = ej if 
fej = kj, and the density pk of infected nodes of degree 
k is given by p k = J^ieK e i = e j > V J e K > where 
N k — P(k)N is the expected number of nodes with degree 
k. Here we have made use of K to denote the set of nodes 



with degree k. This notation allows to group the sums by 
the degrees of the nodes, for instance 

J^ajMkJ 1 )^ = kY,P(k'\k)R x {k'- x )p k , (19) 



(20) 



(21) 



After some algebra eq. ([10]) leads to the generalized HMF 
equation 

= -p Pk + (3k(l - Pk ) J2 Pik'WRxik'-^pv 
k' 

+ /* j Sfcp fc £p(fc / |fc)iZ A (fc , - 1 )ptf 

k' 

+ \p 2 kY J R^k'- 1 ) 2 P{k'\k)p 2 , 

k' 

~\^ k2 (Y,Rx(k'^)P(k'\k) Pk ^ 
If A = 1 then i?i(fc _1 ) = ~ and eq. (J20J becomes 
= + 0k{l - p k ) ^ ^P{k'\k) Pk , 

k' 

+ ^k Pk ^p{k'\k) Pk , + \fk J2 -^p(k'\k)p 2 

k< k' 

If A — > oo then i?oo(fc _1 ) = 1 and eq. (|20| reads 
= -fip k + j3k(l - p k ) ^ P{k'\k) Pk ' 

k' 

+ pPk Pk J2 P(k'\k) Pk ' + l -P 2 kY,P(k'\k)p 2 , 

k> k' 

-^ 2 k 2 (^2P(k'\k)p k ^j (22) 

Again, the first two terms in both cases correspond to the 
standard CP and RP HMF equations respectively, and 
the additional terms arc second order contributions corre- 
sponding to reinfections and multiple infections. 

Conclusions. — We have proposed a new framework 
to study disease spreading in networks. By defining a set 
of discrete-time equations for the probability of individ- 
ual nodes to be infected, we construct a dynamical sys- 
tem that generalizes from an individual contact process 
to the classical case in which all connections are concur- 
rently used, for any complex topology. Solving the equa- 
tions at the stationary state, we find the whole phase dia- 
gram of the system. The numerical solution of the analytic 
equations overcomes the computational cost of MC sim- 
ulations. Moreover, the formalism allows to gain insight 
on the behavior of the critical epidemic threshold for dif- 
ferent values of the probability of contacting a fraction 
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A of neighbors per time step. The proposed model deals 
with infections driven by direct contacts between nodes, 
but not with traffic situations where nodes transmit the 
epidemics by flow communication with others [15] . In this 
latter case, the routing protocol of traffic between nodes 
is absolutely relevant and can change the critical point 
of the epidemic spreading. We are currently working to 
adapt the present formalism also to traffic situations. 
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